function [r_4, moments4, parameters4] = cmd_estimation_prgm(data,ln_y1,ln_y2,ln_a,ln_ri,ln_rc,sa1,ln_a1, ln_a2, ln_r1, ln_r2, m, outname, figout, txtout)
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (1) ENVIRONMENT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Change to user directory
    cd '~\downloads\cmd_replication_public\matlab';   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% (2) COMPUTE SAMPLE MOMENTS %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Vector of sample moments
    mhat_est        = mhat(ln_y1,ln_y2,ln_a,ln_ri,ln_rc);

    % Vector of covariance sample moments from SCF (for 4 parameter model)
    mhat_scf_est    = mhat_scf(ln_r1,ln_r2,ln_a1,ln_a2);
    
    % Vector of additional sample models from SCF (for 4 parameter model)
    mhat_scf_param4_est = mhat_scf_param4(ln_r1,ln_r2,ln_a1,ln_a2);
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% (3) SOLVE CMD MINIMIZATION PROBLEM WITH IPOPT %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    invV = eye(20); % Regardless of spec, we need this

    % Estimate 4-parameter vector
    thetahat4  = cmd_ipopt4(mhat_est,mhat_scf_est,sa1,invV,mhat_scf_param4_est);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% (4) EVALUATE MODEL MOMENTS AT ESTIMATED PARAMETERS %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Vector of model moments evaluated at the estimated 4-parameter vector
    m4         = estmoms(thetahat4,4,mhat_scf_est,sa1,mhat_scf_param4_est);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% (5) BACK OUT PARAMETERS OF INTEREST %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Compute rates of return implied by the 4-parameter model
    [r1_4,r2_4,coef_ri_grp1, coef_rc_grp1, coef_ri_grp2, coef_rc_grp2] = computerates(thetahat4,4,ln_ri,ln_rc);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% (6) COMPUTE STANDARD ERRORS WITH BLOCK BOOTSTRAP %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Data matrix
    data_bs     = [ln_y1,ln_y2,ln_a,ln_ri,ln_rc];

    % SCF data matrix
    scf_bs      = [ln_r1,ln_r2,ln_a1,ln_a2];

    % Length of (overlapping) bootstrap blocks; change from 4 to 3 after
    %   reducing time period covered to 1989-2016
    l = 3;
    
    % Compute standard errors in the 4-parameter model
    [se_thetahat4,se_r1_4,se_r2_4,se_coef_ri_grp1,se_coef_rc_grp1,se_coef_ri_grp2,...
        se_coef_rc_grp2] = stderr(data_bs,scf_bs,4,l,m,sa1,invV, mhat_scf_param4_est);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% (7) PLOT AND EXPORT GROUP-SPECIFIC RATES OF RETURN, %%%%%%%%%%
%%%%%%%%%%%%%%% AS WELL AS POINT ESTIMATES OF PIS, COEFFICIENTS %%%%%%%%%%%
%%%%%%%%%%%%%%% ON INTEREST RATE AND CREDIT RISK TERMS, AND LOGS %%%%%%%%%%
%%%%%%%%%%%%%%% OF INTEREST RATE AND CREDIT RISK TERMS %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Sample size
    n = size(data,1);

    % Compute bounds of 95% confidence intervals in the 4-parameter model
    [lb_r1_4,ub_r1_4] = computeci95(r1_4,se_r1_4,n,4);
    [lb_r2_4,ub_r2_4] = computeci95(r2_4,se_r2_4,n,4);

    % Stack group-specific rates of return and bounds of confidence intervals
    r_4=[data.year,r1_4,r2_4,ub_r1_4,lb_r1_4,ub_r2_4,lb_r2_4,se_r1_4,se_r2_4];
       
    % Plot group rates of return from four-parameter exercise
    figure
    plot(r_4(:,1),r_4(:,2),'Color','blue','Linewidth',1.5)
    hold on
    plot(r_4(:,1),r_4(:,3),'Color','red','Linewidth',1.5)
    hold on
    plot(r_4(:,1),r_4(:,4),'--','Color','blue','Linewidth',1.25)
    hold on
    plot(r_4(:,1),r_4(:,6),'--','Color','red','Linewidth',1.25)
    hold on
    plot(r_4(:,1),r_4(:,5),'--','Color','blue','Linewidth',1.25)
    hold on
    plot(r_4(:,1),r_4(:,7),'--','Color','red','Linewidth',1.25)
    hold on
    xlabel('Year','interpreter','latex')
    ylabel('Rate of Return','interpreter','latex')
    legend('Group 1','Group 2',...
    'Location','northeast',...
    'interpreter','latex','NumColumns',1)
    legend boxoff
    set(gca,'TickLabelInterpreter','latex')
    print(fullfile(figout,['plot-rates4-btstrp',sprintf('%d',m),'-',outname]),...
    '-djpeg','-r500')
    hold off

    % Export data
    newdirname = strcat(outname, '-btstrp', sprintf('%d',m));
    newdir = fullfile(txtout,newdirname);
    mkdir(newdir);
    
    csvwrite(fullfile(newdir,'r4.txt'), r_4);
    
    moments4 = [mhat_est, m4];
    csvwrite(fullfile(newdir,'fitted_moments4.txt'),moments4);
    
    parameters4=[thetahat4, se_thetahat4];
    csvwrite(fullfile(newdir,'parameters4.txt'),parameters4);
    
    coef_ses = [coef_ri_grp1, se_coef_ri_grp1; coef_rc_grp1, se_coef_rc_grp1; ...
                coef_ri_grp2, se_coef_ri_grp2; coef_rc_grp2, se_coef_rc_grp2];
    csvwrite(fullfile(newdir,'coef_ses4.txt'),coef_ses);
    
    scfcalibrated4 = [mhat_scf_param4_est; sa1];
    csvwrite(fullfile(newdir,'scfcalibrated4.txt'),scfcalibrated4);
end